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ABSTRACT 

A simple dynamical model is employed to study the possible orbital evolution 
of scattered planets and phase plane analysis is used to classify the parameter 
space and solutions. Our results reconfirm that there is always an increase in 
eccentricity when the planet was scattered to migrate outward when the initial 
eccentricity is zero. Applying our study on the Solar System and considering 
the existence of the Kuiper Belt, this conclusion implies that Neptune was 
dynamically coupled with the Kuiper Belt in the early phase of the Solar 
System, which is consistent with the simulational model in Thommes, Duncan 
& Levison (1999). 

Subject headings: celestial mechanics - stellar dynamics - planetary systems 
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1. Introduction 

Due to the observational effort, both the number of extra-solar planets and small 
bodies in the outer Solar System has increased dramatically. These new discoveries provide 
interesting laboratories to study theories of planetary dynamics (Jiang & Ip 2001). On the 
other hand, some of the observational phenomena require explanation from the dynamical 
models. For example, in order to explain the dynamics of Plutinos, the possible migration 
of Neptune's orbit has been suggested by Malhotra (1995). Fernandez & Ip (1984) proposed 
that the interaction with planetesimals might make the orbit of Neptune expand and 
therefore migrate slowly. This "slow migration" was demonstrated by Hahn & Malhotra 
(1999). 

On the other hand, Thommes, Duncan & Levison (1999) successfully showed that the 
Neptune might form in the Jupiter-Saturn region of the Solar System and was scattered to 
migrate outward to the current location due to the interaction with the growing Jupiter or 
Saturn. In their model, Neptune was scattered to have a "fast migration" at first and then 
transfer to "slow migration" due to the interaction with planetesimals. 

Numerical simulations provide opportunities for us to explore the possible solutions for 
astronomical problems and are always good tools for more practical modelling. Nevertheless, 
it is important to use the analytical model to classify the solutions on the parameter space 
and to understand the principle physical processes within the simulations. 

Considering two planets moving around a central star, if the apocentre-distance of the 
inner planet is much smaller than the pericentre-distance of the outer planet, there would 
be no scattering between these two planets. They can only influence each other through 
secular perturbation. Therefore, the scattering might happen when the inner planet (say, 
the smaller one) moving on a more eccentric orbit and the outer planet (say, the giant 
planet) moving on a circular orbit and their orbits cross each other. We can call this process 
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"crossing scattering" and the Tisserand relation is valid if the distance between these two 
planets never gets too close. 

Another type of scattering process might happen for two planets both moving on 
circular orbits, as the newly formed Neptune and Saturn in Thommes, Duncan & Levison 
(1999). Let us assume the small planet's orbital radius (r s ) is slightly shorter than the 
giant planet's orbital radius (r g ) and so the angular velocity of the small planet (dO s /dt) is 
slightly larger than the giant planet's angular velocity (dO g /dt). If both of them moving 
anti-clockwise as the convention of a polar coordinate, there must be a chance that they 
approach to each other and go through three periods: (1) $ s — 9 g — e, (2) 9 S = 9 g and (3) 
9.s = 9 g + e, where e is a small number. During the first period, the small planet should 
gain the angular momentum from the giant planet because the giant planet is pulling 
the small planet in the direction which the small planet is moving to. Similarly, during 
the final period, the small planet should lose the angular momentum to the giant planet 
because the giant planet is pulling the small planet in the direction oppositive to the small 
planet's moving direction. Obviously, during the middle period, there should be no angular 
momentum exchange and the interaction between two planets is stronger because they are 
closer to each other. The net effect of this type of scattering is that the small planet is 
pulling away from the central star and the angular momentum is approximately conserved. 
We can call this process "swing scattering" . 

If one assume that all planets are moving on circular orbits initially as in Thommes, 
Duncan & Levison (1999), this "swing scattering" definitely happen before the ordinary 
"crossing scattering" could take place. Therefore, both of them are important for the overall 
orbital evolution. 

In this paper, we will focus on "swing scattering" only. We analytically model this 
process by adding a forcing term on the gravitational force from the central star. This term 
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represents the interaction with the growing giant planets. We assume that the scattering is 
in the radial direction, so the angular momentum is still conserved for the planet. 

We study the possible outcome under this scattering. Thus, in our system, the planet 
would face three stages: (1) pre-scattering stage: freely orbiting a central star in a two-body 
dynamical system. (2) scattering stage: influenced by the scattering term in additional to 
the force from the central star. (3) post-scattering stage: the scattering stops and return to 
a two-body system again. 

In Section 2, we describe the equations for our model. We do the phase plane analysis 
in Section 3 and report the results in Section 4. In Section 5, we provide the conclusions. 

2. The Model 

We consider a two-body system that a planet is moving around a central star on 
the two dimensional space governed by the attractive inverse square law of force with an 
additional term representing the scattering. Therefore, 

/(r) = ^ + *i, (1) 

where k is positive and the minus sign ensures that the force is toward the central star. 
Moreover, the parameter k\ is also positive, representing the scattering force in the direction 
away from the star. We assume that k\ is a constant during the scattering. 

Let u = 1/r, from Goldstein (1980), we have the following equation 

<P9 ~ u Z 2 u 2 ' K) 

where m is the mass and / is the angular momentum of the planet. We use the polar 
coordinate (r, 9) to describe the location of the planet. 
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The angular momentum I is a constant here, so we have 

mr 2 d6 = Idt. (3) 

Because of this, we can use 9 as our independent variable. We use 6 to label time t afterward 
and one can easily gets t from the above equation. 

From Equation (|1|) and Equation (0), we have 

d 2 u _ $ 



where /3 = mk/(l 2 ) and /3i = mki/(l 2 ). One can see that u, /3, and /?i are all positive. 

We then transform Equation (^) to the following autonomous system: 

du _ 

d9 = V = ^M 



(5) 



— = -u + /3 - = f 2 (u,v) 

du u z 



In order to understand how the solutions depend on the parameters, we will do the 
phase plane analysis to classify the solutions completely in terms of (3 and fli. 

During the pre-scattering and post-scattering stage when (3\ = 0, there is one fixed 
point (u, v) = (j3, 0) on the phase plane. The eigenvalues corresponding to this fixed point 
are complex, i.e. A = ±i, so it is a stable center point. In the next section, we begin to do 
the analysis for the scattering stage when /3i > 0. 

3. Phase Plane Analysis 

The fixed points (u, v) of Equation (|5|) satisfy the following equations: 

v = Q, 

u 3 - (3u 2 + fa = 0, (6) 
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where «^0. Let 



equation (|6|) becomes 
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The solutions of the above cubic equation completely depend on A (Neumark 1985). 
We will consider all different values of A except there are a few cases which are not inside 
our scope. We will not consider the case when A < — 1 because this gives us f3\ < but we 
assume (3\ > in this paper. On the other hand, the probability that A = 1 and A = are 
very small because the scattering parameter (3\ need to precisely equal to particular number 
for a given (3. We will not discuss the case of A = 1 because the eigenvalues are zero for 
the first order linearization of the phase plane analysis. The results for A = will be in 
Appendix A. 

In fact, the value of A has the physical meaning of the strength of the scattering 
relative to the strength of the central star's gravity. Thus, A = — 1 means there is no 
scattering and the larger A corresponds to the stronger scattering relatively to the gravity 
of the central star. 

Therefore, in this section, there would be three cases, which correspond to different 
relative strength of the scattering : A > 1 (27ft - A(3 3 > 0), < A < 1 (2/3 3 < 27ft < 4/3 3 ) 
and -1 < A < (0 < 27ft < 2/3 3 ). 



3.1. Case 1: A > 1 

In this case, the solution of equation (§) is 

y = cosh?/, (10) 

where r\ satisfies 

27/ 
cnsh 3?i = A = 

2/5 



cosh3r / = A= 27/Jl ~ 2/j3 . (11) 



Equation (0) gives 

71 = 

3 
Therefore, the fixed point is 



u = ^(-2y+l). (12) 



(u,v)= f|(-2cosh7; + l),0J. (13) 

It is easy to show that u is negative for this fixed point. The solution curves in the 
neighborhood of this fixed point will be unphysical. 

3.2. Case 2: < A < 1 

For this case, the solutions of equation (||) are 

yi = cos f, (14) 

and 

2/2,3 = - cos (- ± n , (15) 



where £ satisfies 

27ft- 2/3 3 
"2^ 

and < £ < |. Therefore, by Equation (|i~2f) , the three fixed points are 



cos3£ = A = ^— (16) 



(«i,«i)= (|(-2cos£ + l),0j (17) 



and 
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(u 2 ,v 2 )= ^(2cos(| + + l),o], (18) 



(u 3 ,vs)= ^(2cos(|-0 + l),o]. (19) 



It can be shown that U\ is negative, so the solution curves near the first fixed point, 
(■Ui,t>i), is unphysical. 

We can also prove that the second fixed point, (^2,^2), is an unstable saddle point and 
the third fixed point is a stable center point (See Appendix B). 

3.3. Case 3: -1 < A < 

The solutions of equation (|8|) are 

2/i = -cos0, (20) 

and 

7/r, Q = nns I 



2/2,3 = COS (^±0), (21) 



where <p satisfies 

27 / 3j_- 2/? 3 

and < 4> < f . Therefore, the fixed points are 



cos30=— A = — r (22) 



and 



(«i,vi) = (^(2cos0 + l),O), (23) 

(«2,v 2 ) = (f(-2cos(| + 0) + l),O), (24) 



(us, v 3 ) = (|(-2 cos(^ - 0) + 1), 0). (25) 
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We can prove that the first fixed point, (ui,Vi), is an stable center point and the second 
fixed point, (u2,v 2 ) is an unstable saddle point (See Appendix C). 

It can be easily shown that U3 is negative, so the solution curves near the third fixed 
point, (-U3,f 3 ), is unphysical. 

4. The Results 

The phase plane analysis in the last section classifies the solution curves by the relative 
strength parameter A (or parameters /3 and Pi). 

There are therefore plenty of possible orbits during the scattering, especially for Case 
2 and Case 3. Because both Case 2 and 3 have one stable center and one unstable saddle 
point in the u > region, we only need to use the results in Case 2 to demonstrate the 
possible migrations of the planet in this section. 

Fig.l shows that a planet initially moving on a orbit with semi-major axis a = 0.94 
can migrate to a = 1.57 and the eccentricity increases from e = 0.33 to e = 0.68. (One can 
regard the unit of length to be AU in this paper.) The u — v plane in Fig. 1(a) shows that 
there are two fixed points: one is an unstable saddle point and another is a stable center 
point. Fig.l(b)-(d) are the orbital evolution on the x — y plane and the time variation 
for both the eccentricity and semi-major axis for the planet. The solid line is for the 
pre-scattering stage, the dotted line is for the scattering stage and the dashed line is for the 
post-scattering stage. The scattering comes in when 9 = 2n and leaves when 9 = 3tt, so 
the scattering acts for half of the orbital period. The dotted line in Figl(a) is the solution 
curve used for the scattering stage in Fig.l(b)-(d). 

In order to simulate the scattering which takes place at larger semi-major axis, we 
choose the value of P to make the initial semi-major axis of the planet to be around 24 
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and the initial eccentricity to be zero. Fig. 2 shows that the planet initially moving on a 
orbit with semi-major axis a = 23.8 can migrate to a = 24.8 and the eccentricity increases 
from e = 0.0 to e = 0.2. Fig.2(b)-(d) are the orbital evolution on the x — y plane and the 
time variation for both the eccentricity and semi-major axis for the planet. The scattering 
begins when 9 = 2n and stops when 9 = 3ir, so the scattering acts for half of the orbital 
period. The dotted line in Fig. 2 (a) is the solution curve used for the scattering stage in 
Fig. 2 (b)-(d). The planet does not migrate as much as in Fig.l because the solution curve 
during the scattering stage is much closer to the fixed point and the variation in u (and 
therefore r) is much smaller. 

To model a migration caused by a stronger scattering, we now assume that there is 
velocity discontinuity when the planet enters the scattering stage from the pre-scattering 
stage. Fig. 3 shows that a planet initially moving on a orbit exactly the same as the initial 
orbit in Fig.2 can migrate to a = 69.0 and the eccentricity increases from e = 0.0 to 
e = 0.81. Fig.3 (b)-(d) are the orbital evolution on the x — y plane and the time variation 
for both the eccentricity and semi-major axis for the planet. 

The scattering starts when 9 = 2n and stops when 9 = 2.55ir so the scattering acts 
for about quarter of the orbital period. The dotted line in Fig3(a), the u — v plane, is the 
solution curve used for the scattering stage in Fig.3(b)-(d). 

We found that the eccentricity always increases and does not decay to zero. In this 
model, the only possibility to make the planet move on a orbit with zero eccentricity again 
is if the scattering stage stops precisely when the solution curve arrives at the fixed point of 
the post-scattering stage. This case can be neglect because the probability is too small. 

Because in the simulational model of Thommes, Duncan and Levison (1999) the 
scattering takes place for many times before Neptune leaves the Jupiter-Saturn region, one 
might want to generalize our model to multiple scattering. One can choose small f3\ for each 
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process and set /3i — between two scattering. All the results of phase plane analysis are 
still the same and one only needs to pick up the correct solution curves to use and connect 
them together. 

5. Conclusions 

We have studied the dynamics of a planet under the influence of "swing scattering" . 
The governing equations can be transformed to an autonomous system and the standard 
phase plane analysis was used to classify the parameter space and possible solutions. 

When A > 1, i.e. f3\ > 4[3 3 /27, there is only one fixed point, where u is negative. The 
solution curves near this fixed point is therefore unphysical. 

When < A < 1, i.e. 2/3 3 /27 < j3i < 4/3 3 /27, there are two fixed points where u > 0. 
The two eigenvalues corresponding to u 2 are real with opposite sign, so u 2 is an unstable 
saddle point. The two eigenvalues corresponding to Uj, are pure complex with opposite sign, 
so -u 3 is a stable center point. 

When — 1 < A < 0, i.e. < j3± < 2(3 3 /27, there are two fixed points where u > 0. The 
two eigenvalues corresponding to U\ are pure complex with opposite sign, so u\ is a stable 
center point. The two eigenvalues corresponding to u 2 are real with opposite sign, so u 2 is 
an unstable saddle point. 

Thus, we found that the solution topology is the same when — 1 < A < 1 (with stable 
fixed point at physical region) but becomes very different when A > 1 (without stable fixed 
point at physical region) 

Physically, this means when A > 1, the scattering is so strong such that the planet 
will be forced to leave this system if the scattering continues, no matter what position 
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on the phase plane is when the scattering happens. When — 1 < A < 1, the scattering is 
weaker and the planet has chances to stay in this system even if the scattering continues, 
depending on what position on the phase plane is when the scattering begins. 

Therefore, A > 1 is an important instability criteria for "swing scattering" . 

These analytical results enable us to study all possible orbits of outward planetary 
migration due to "swing scattering". We found that the planet can migrate from an orbit 
with non-zero eccentricity to another eccentric orbit for both inner and outer part of the 
planetary system. In some situation, when the scattering is stronger, the planet (or small 
bodies like Kuiper Belt Objects) can migrate to 100 AU. 

Our results reconfirm that the migration caused by scattering always involves the 
eccentricity increasing during the orbital evolution. There is only one possibility that the 
orbit can be completely circular again after the scattering stage. That is if the scattering 
stage stops precisely when the solution curve arrives at the fixed point of the post-scattering 
stage. We can neglect this because the probability is extremely small. 

On the other hand, the orbit of the planet (Neptune) could be circularized if it interact 
with the disc (Kuiper Belt) in the outer planetary system (Solar System) after the orbital 
scattered migration as we see in the simulations of Thommes, Duncan & Levison (1999). 

If this was the process for the outer Solar System, i.e. Neptune was scattered to 
migrate outward but increased the eccentricity at first, then circularized by the Kuiper 
Belt afterward, the Kuiper Belt has to be massive enough. This might indirectly supports 
the picture that the Kuiper Belt was at least two order more massive: an argument from 
the collision and formation history of the Kuiper Belt. If that is true, a large fraction 
of the Kuiper Belt Objects (KBOs) would escape from the region of 30-50 AU and part 
of those non-escaped KBOs would be captured by the Neptune into 3:2 resonance. The 
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more complete physical picture might be that the Neptune and the massive KBOs were 
coupled dynamically. They influenced each other, i.e. KBOs circularized the Neptune's 
orbit and the Neptune scattered some of the KBOs, until most KBOs left and the remained 
KBOs can then be treated as test particles under the influence of the outward migration 
of Neptune. This picture might be important for further study of the orbital evolution for 
both Neptune and KBOs. 
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A. The Case when A = 



From (RJ), we know the eigenvalues satisfy 

(§H (§H - (f ) (f ) - 

Since fi(u, v ) = v and /2(w, f) = — u + j3 — @k, the above equation becomes 

A 2 - (-1 + ^)=0. (A2) 



Thus 



u 3 



* - 4 1 + C' 
- ±(^*f ■ 

The solutions of Equation (||) are 

yi = 0, (A4) 

and 



1/2,3 = ±^- (A5) 
Therefore, from Equation (|i2|), three fixed points are 

(«i,«i) = (f,0), (A6) 

(«a,«Si) = (|(l + V3),0), (A7) 

(«8,«s) = (^(l-V3),0). (A8) 

From Equation (|A3|) , the eigenvalues corresponding to («i,i>i) are ±V3 and the 



eigenvalues corresponding to (^2,^2) are ±J—1 + . /-. 3 . Thus, fixed point (ui,vx) is an 
unstable saddle point, and (^2^2) is a stable center point. 

Since W3 is negative, the solution curve in the neighborhood of (1*3,^3) is unphysical. 
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B. The Eigenvalues for Case 2 

From Equation (|T6|), we have 



f3 1 = l-P* C o S 3Z + ^f3 3 . (Bi; 



Let 



H x &0) = u 2 -(2^y 



^(2cos(|+0 + l)-f(4cos3^ + 4)3 
|{ 2o os(| + + l-(4cos(3fl+4)i 



Thus, 



dH x f3 



Numerically, we can easily show that 



|{-2sin(| + 0H-4sin3£(4cos3£ + 4)-i}. (B2) 



- 2 sin(- + 0+4 sin 3£(4 cos 3£ + 4)~§ < 0, (B3) 

o 

for < £ < |. Further, because (3 > 0, we then have 

We also have #i(0, /?) = 0, so it implies #i(£, /3) < 0. Therefore, u| - 2ft < for all /? > 
and < £ < | . 

Thus, from Equation (|^), the two eigenvalues corresponding to «2 are real with 
opposite sign, so it is an unstable saddle point. 

Since < f < § and (&)s < | < (f )s, we have 



2 <2cos(^-f) + l < V3 + 1. 



-17- 

Thus, 

u 3 = I (2cos(^ - + l) > ~/3 > 2(^)i = (2ft)*. 

Thus, W3 — 2/?! > 0. From Equation ([A 3D, the two eigenvalues corresponding to u 3 are pure 
complex with opposite sign, so it is a stable center point. 

C. The Eigenvalues for Case 3 

Since < ft < ^/? 3 and < < | in this case, it implies 

3 3 
w?-2ft = ^(2cos0+l) 3 -2ft 

> ^{(2cos0+l) 3 -4} 

> ^{(Vs + lf-4} 

> 



Therefore, from Equation (|A3|), the two eigenvalues corresponding to U\ are pure 



complex with opposite sign, so it is a stable center point. 
From Equation fl2"2]), we have 



/3i = -^ 3 cos30+^/3 3 . (CI) 



Let 



tf 2 (0,ft = w 3 -2/? ; 



:■! 



= (|J (_2cos(^ + 0) + l) 3 -2ft 

= ^-{(-2cos(^ + 0) + l) 3 + 4cos(30)-4}. 

Follow the same approach as in Appendix B, we can show that ^(0, /3) < 0, for all 
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(3 > and < < |. Therefore, we have u\ — 2(5\ < 0. From Equation (|A3|), the two 
eigenvalues corresponding to u 2 are real with opposite sign, so it is an unstable saddle point. 
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Fig. 1.— The Result when /3 = 1.2, f3\ = 2.2(3 3 /27. (a) The u — v phase plane, where the 
dotted line is the solution curve used during the scattering stage, (b) The orbit on the x — y 
plane, (c) The eccentricity as function of 9, (d) The semi-major axis as function of 9, where 
the solid line is for the pre-scattering stage, the dotted line is for the scattering stage and 
the dashed line is for the post-scattering stage. 
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Fig. 2.— The Result when f3 = 0.04, f3\ = 2.2/5 3 /27. (a) The u — v phase plane, where the 
dotted line is the solution curve used during the scattering stage, (b) The orbit on the x — y 
plane, (c) The eccentricity as function of 9, (d) The semi-major axis as function of 9, where 
the solid line is for the pre-scattering stage, the dotted line is for the scattering stage and 
the dashed line is for the post-scattering stage. 
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Fig. 3.— The Result when f3 = 0.04, f3\ = 2.2/5 3 /27. (a) The u — v phase plane, where the 
dotted line is the solution curve used during the scattering stage, (b) The orbit on the x — y 
plane, (c) The eccentricity as function of 9, (d) The semi-major axis as function of 9, where 
the solid line is for the pre-scattering stage, the dotted line is for the scattering stage and 
the dashed line is for the post-scattering stage. 



